#include "stdafx.h"

#include "Solver.h"
#include <math.h>
#include "matrix.hpp"
#include "evector.hpp"
Solver::Solver()
{
    //ctor
}

Solver::~Solver()
{
    //dtor
}


EVector<double> Solver::CGMethod(FullMatrix<double> B, EVector<double> F, double tol)
{
	int degree = F.size();
    EVector<double> _x(degree);
    EVector<double> _r(degree);
    EVector<double> _p(degree);
    EVector<double> _w(degree);
    EVector<double> temp(degree);
	FullMatrix <double> Bp(degree,degree);
	_x=F;
	_r=F-B*_x;
	_p=_r;
	double rho0,rho1,anpha;
	rho0=(_r *_r).GetNorm();
	for(int i=1;i<=100;i++)
	{
	_w = B*_p;
	anpha=rho0/(_p*_w).GetNorm();
	_x=_x+anpha*_p;
	_r=_r-anpha*_w;
	rho1=(_r*_r).GetNorm();
	if(sqrt(rho1)<tol)
		break;
	else
		_p=_r+rho1/rho0*_p;
	rho0=rho1;

	
	}
	return _x;
  
}
EVector<double>Solver:: giaiHePT(FullMatrix<double> B, EVector<double> F)
{
	int n = B.getNrCols();
	EVector<double> tmp(n);
	tmp.setZero();
	for (int i = 1; i<n-1; i++)
	{
		B[i-1][i] /= B[i-1][i-1];
		F[i-1] /= B[i-1][i-1];
		B[i-1][i-1] = 1.0;
		for (int j = i; j<=i+1;j++)
		{
			B[i][j] -= B[i][i-1]*B[i-1][j];
		}
		F[i] -= B[i][i-1]*F[i-1];
		B[i][i-1] = 0;
	}
	if (n>=2)
	{
		B[n-2][n-1] /= B[n-2][n-2];
		F[n-2] = F[n-2]/ B[n-2][n-2];
		B[n-2][n-2] = 1.0;
		B[n-1][n-1] -= B[n-2][n-1]*B[n-1][n-2];
		F[n-1] -= B[n-1][n-2]*F[n-2];
		B[n-1][n-2] = 0;
	}
	F[n-1] /= B[n-1][n-1]; B[n-1][n-1] = 1.0;

	for (int i = n-2; i>=0; i--)
	{
		F[i] -= F[i+1] * B[i][i+1];
		B[i][i+1] = 0;
	}
	return F;
}